In the present report, we first estimate the narrow-sense heritability and evolvability at the species level, i.e. across all populations. Then we estimate these two parameters at the population level, i.e. one parameter value per population. For that, we first run one model per population (hereafter referred as ‘population-specific models’). Second, we run a full model with all populations (hereafter referred as ‘full model’), based on Archambeau et al. (2023).
As heritability and evolvability are specific to the environment, the models are fitted in each common garden separately. We fit models that do not account for the nursery effect for the two common gardens in which trees come from the same nursery: the sites of Yair (FS) in the Scottish borders (in which all trees come from NG) and Glensaugh (FE) in which all trees come from NE (except four trees that come from NG and that we remove for the analyses).
We fit models that do account for the nursery effect in the common garden in which trees come from different nurseries: Inverewe (FW), which contains cohorts of trees raised in each of the three nurseries as follows: 290 grown locally in the NW; 132 grown in the NG; and 82 grown in the NE.
The data
For details on the data, see report EstimatingPhenotypicPlasticityScotsPine.qmd.
Code
# We load the field datadata <-read_excel(here("data/Field.xlsx"), na ="NA")# we load the nursery data (dataset updated by Annika - 02/06/2023)nursery_data <-read_excel(here("data/SPnurseries_updatedByAnnika02062023.xlsx")) %>% dplyr::rename(FieldCode=Tag)# we merge the twodata <- data %>%left_join(nursery_data, by=c("FieldCode"))
We have to change the names of the block to specify that they are nested within sites:
Code
# site-specific blocksdata <- data %>%mutate(Block =paste0(FieldSite,"_",Block))
According to Beaton et al. (2022), there are four randomised blocks in both FS and FE and three in FW. We look at the number of individuals in each block:
Code
data %>%group_by(Block) %>%summarise("Number of individuals"=n()) %>%kable_mydf()
Block
Number of individuals
FE_A
168
FE_B
168
FE_C
168
FE_D
168
FS_A
168
FS_B
168
FS_C
168
FS_D
168
FW_A
168
FW_B
168
FW_c
1
FW_C
167
We see that there is one typo in FW (block noted as c instead of C), that we have to correct:
Code
data <- data %>%mutate(Block=case_when(Block=="FW_c"~"FW_C",TRUE~ Block))data %>%group_by(Block) %>%summarise("Number of individuals"=n()) %>%kable_mydf()
Block
Number of individuals
FE_A
168
FE_B
168
FE_C
168
FE_D
168
FS_A
168
FS_B
168
FS_C
168
FS_D
168
FW_A
168
FW_B
168
FW_C
168
In the present report, the estimation of heritability and evolvability is done for:
Code
trait <-"HA20"
We subset the dataset by removing missing values for the trait considered and by removing the four trees in FE that belong to the NG nursery:
Code
site_codes <-unique(data$FieldSite)subdata <- data %>% dplyr::select(FieldSite,Block,PopulationCode,Family,Nursery,all_of(trait)) %>%drop_na() %>% dplyr::rename(site=FieldSite,bloc=Block,pop=PopulationCode,fam=Family,nurs=Nursery,y=any_of(trait)) %>%filter(!(site=="FE"& nurs =="NG")) # we remove the four trees that come from the nursery NG in the field site FE
Species-level parameters
In this section, we estimate evolvability (\(I\)) and narrow-sense heritability (\(h^2\)) at the species level, i.e. one value for all populations.
Model equation
Here is the mathematical model of the model for a common garden with different nurseries. The mathematical model of the model for a common garden with trees coming from one nursery is the same, without the nursery intercepts.
\(\beta_{0}\) is the global intercept. \(B_{b}\) are the block intercepts with a variance \(\sigma^{2}_{B}\). \(P_{p}\) are the population intercepts with a variance \(\sigma^{2}_{P}\). \(F_{f}\) are the family intercepts with a variance \(\sigma^{2}_{F}\). \(N_{n}\) are the nursery intercepts with a variance \(\sigma^{2}_{N}\). \(\sigma^{2}_{r}\) is the residual variance.
Partitioning of the total variance \(\sigma_{tot}^{2}\):
S4 class stanmodel 'anon_model' coded as follows:
// Variance partitioning in the common gardens with trees from different nurseries
data {
int<lower=1> n; // Number of observations
vector[n] y; // Response variable
int<lower=0> n_bloc; // Number of blocks
int<lower=0> n_pop; // Number of populations
int<lower=0> n_fam; // Number of families
int<lower=0> n_nurs; // Number of nurseries
int<lower=0, upper=n_bloc> bloc[n]; // Blocks
int<lower=0, upper=n_pop> pop[n]; // Populations
int<lower=0, upper=n_fam> fam[n]; // Families
int<lower=0, upper=n_nurs> nurs[n]; // Nurseries
}
parameters {
real beta0; // global intercept
simplex[5] pi;
real<lower = 0> sigma_tot;
vector[n_bloc] z_block;
vector[n_pop] z_pop;
vector[n_fam] z_fam;
vector[n_nurs] z_nurs;
}
transformed parameters {
real R_squared;
real<lower = 0> sigma_r;
real<lower = 0> sigma_block;
real<lower = 0> sigma_pop;
real<lower = 0> sigma_fam;
real<lower = 0> sigma_nurs;
vector[n_bloc] alpha_block;
vector[n_pop] alpha_pop;
vector[n_fam] alpha_fam;
vector[n_nurs] alpha_nurs;
vector[n] mu; // linear predictor
// variance partitioning with the simplex pi
sigma_r = sqrt(pi[1]) * sigma_tot;
sigma_block = sqrt(pi[2]) * sigma_tot;
sigma_pop = sqrt(pi[3]) * sigma_tot;
sigma_fam= sqrt(pi[4]) * sigma_tot;
sigma_nurs= sqrt(pi[5]) * sigma_tot;
alpha_block = z_block*sigma_block;
alpha_pop = z_pop*sigma_pop;
alpha_fam = z_fam*sigma_fam;
alpha_nurs = z_nurs*sigma_nurs;
mu = rep_vector(beta0, n) + alpha_block[bloc] + alpha_pop[pop] + alpha_fam[fam] + alpha_nurs[nurs];
R_squared = 1 - variance(y - mu) / variance(y);
}
model{
//Priors
beta0 ~ normal(mean(y),20);
sigma_tot ~ student_t(3, 0.0, 1.0);
z_block ~ std_normal();
z_pop ~ std_normal();
z_fam ~ std_normal();
z_nurs ~ std_normal();
// Likelihood
y ~ normal(mu, sigma_r);
}
generated quantities {
//Variances, narrow-sense heritability and evolvability
real<lower=0> sigma2_tot;
real<lower=0> sigma2_r;
real<lower=0> sigma2_block;
real<lower=0> sigma2_pop;
real<lower=0> sigma2_fam;
real<lower=0> sigma2_nurs;
real<lower=0> h2;
real<lower=0> I;
sigma2_tot = square(sigma_tot);
sigma2_r = square(sigma_r);
sigma2_block = square(sigma_block);
sigma2_pop = square(sigma_pop);
sigma2_fam = square(sigma_fam);
sigma2_nurs = square(sigma_nurs);
h2 = sigma2_fam/(sigma2_r+sigma2_fam);
I = sigma2_fam/(mean(y)^2);
}
Running the models
Options for the Stan models:
Code
# Sampling in Bayesian modelsn_chains <-4# number of chains (MCMC)n_iter <-2500# number of iterationsn_warm <-1250# number of iterations in the warm-up phasen_thin <-1# thinning intervalsave_warmup <-FALSE# Credible intervalsconf_level <-0.95
Warning: There were 3 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess
Interval plots
Variance partitioning
Code
readRDS(file=here(paste0("outputs/Evolvability/",trait,"_SpeciesLevelModels.rds"))) %>%lapply(function(x){ broom.mixed::tidyMCMC(x,pars=("pi"),droppars =NULL, ess = F, rhat = F, conf.int = T, conf.level = conf_level) }) %>%bind_rows(.id ="site") %>%mutate(prop_var =case_when(term =="pi[1]"~"Residuals", term =="pi[2]"~"Blocks", term =="pi[3]"~"Populations", term =="pi[4]"~"Families", term =="pi[5]"~"Nurseries")) %>%ggplot(aes(y = prop_var, x = estimate, xmin = conf.low, xmax = conf.high, color=site)) +geom_pointinterval(position =position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +ylab("") +xlab("Proportion of variance explained") +labs(color ="Field sites") +theme(axis.text =element_text(size=18),panel.grid.minor.y=element_blank(),panel.grid.major.y=element_blank()) +guides(color=guide_legend(ncol=1))
In this section, we estimateevolvability (\(I\)) and narrow-sense heritability (\(h^2\)) at the population level, i.e. one value for each population.
Full dataset
We first fit the models on the entire dataset.
Population-specific models
We first fit one model per population.
Here is the model equation for a common garden with different nurseries. The model equation for a common garden with trees coming from one nursery is the same, without the nursery intercepts.
\(\beta_{0}\) is the global intercept. \(B_{b}\) are the block intercepts with a variance \(\sigma^{2}_{B}\). \(F_{f}\) are the family intercepts with a variance \(\sigma^{2}_{F}\). \(N_{n}\) are the nursery intercepts with a variance \(\sigma^{2}_{N}\). \(\sigma^{2}_{r}\) is the residual variance.
Partitioning of the total variance \(\sigma_{tot}^{2}\):
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 10 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 7 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 2 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 3 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
\(\beta_{0}\) is the global intercept. \(B_{b}\) are the block intercepts with a variance \(\sigma^{2}_{B}\). \(P_{p}\) are the population intercepts with a variance \(\sigma^{2}_{P}\). \(F_{f}\) are the family intercepts with a variance \(\sigma^{2}_{F}\). \(N_{n}\) are the nursery intercepts with a variance \(\sigma^{2}_{N}\). \(\sigma^{2}_{r}\) is the residual variance.
The prior of \(\beta_{0}\) was weakly informative and centered around the mean of the observed values for the trait under considered, as follows:
\[\beta_{0} \sim \mathcal{N}(\mu_{y},20)\]
The population, nursery and block intercepts, \(P_{p}\), \(N_{n}\) and \(B_{b}\) were considered normally-distributed with variances \(\sigma^{2}_{P}\), \(\sigma^{2}_{N}\) and \(\sigma^{2}_{B}\), such as:
The family intercepts \(F_{f(p)}\) were considered to follow some population-specific normal distributions, such as: \[F_{f(p)} \sim \mathcal{N}(0,\sigma^{2}_{F_{p}})\]
where \(\sigma^{2}_{F_{p}}\) are the population-specific variances among families.
To partition the total variance, we parameterize our model so that only the total variance, \(\sigma_{tot}^{2}\) has a prior, such that:
where \(\overline{\sigma_{F_{p}}}\) and \(\overline{\sigma_{F_{p}}}^{2}\) are the mean of the population-specific among-family standard deviations (\(\sigma_{F_{p}}\)) and variances (\(\sigma^{2}_{F_{p}}\)), respectively, and \(\sum_{l}^{5}\pi_{l}=1\) (using the simplex function in Stan).
The population-specific among-families standard deviations \(\sigma_{F_{p}}\) follow a log-normal distribution with mean \(\overline{\sigma_{F_{p}}}\) and variance \(\sigma^{2}_{K}\), such as:
Warning: There were 68 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: There were 72 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess
Warning: There were 8 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
The models have divergent transitions. We can visualize them with pairs plots.
Code
readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel.rds"))) %>%lapply(function(x){mcmc_pairs(as.array(x), np =nuts_params(x), pars =c("pi[1]","pi[2]","pi[3]","pi[4]"),off_diag_args =list(size =0.75)) })
$FE
$FW
$FS
Comparing \(h^2\) and \(I\) estimates
We check the correlations among \(h^2\) and \(I\) values estimated calculated with population-specific models or with the full model including all populations.
Code
list_df <-lapply(c("I","h2"), function(x){# Extract the median estimates of h2/I in the population-specific models# ======================================================================df_pop_specific_models <-readRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels.rds"))) %>%lapply(function(site_i){ site_i %>%lapply(function(pop_i){ broom.mixed::tidyMCMC(pop_i,pars=(x),droppars =NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) }) %>%bind_rows(.id ="pop") }) %>%bind_rows(.id ="site") %>%mutate(model="pop_specific_models") %>%mutate(model_site=paste0(model,"_",site)) %>% dplyr::select(-term)# Extract the median estimates of h2/I in the population-specific models# ======================================================================df <-readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel.rds"))) %>%lapply(function(site_i){ broom.mixed::tidyMCMC(site_i,pars=(paste0(x,"_pop")),droppars =NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) }) %>%bind_rows(.id ="site") %>%mutate(model="full_model") %>%mutate(model_site=paste0(model,"_",site)) %>%mutate(pop=rep(unique(subdata$pop) %>% sort,3)) %>% dplyr::select(-term) %>%bind_rows(df_pop_specific_models)}) %>%setNames(c("I","h2"))list_df[["h2"]] %>%ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site,shape=model)) +geom_pointinterval(position =position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +xlab("") +ylab("Narrow-sense heritability") +theme(axis.text =element_text(size=18),panel.grid.minor.x=element_blank(),panel.grid.major.x=element_blank()) +guides(color=guide_legend(ncol=1))
Warning: There were 99 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Code
readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds"))) %>%lapply(function(x){mcmc_pairs(as.array(x), np =nuts_params(x), pars =c("pi[1]","pi[2]","pi[3]","pi[4]"),off_diag_args =list(size =0.75)) })
$FE
$FS
Code
readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds"))) %>%lapply(function(x){mcmc_pairs(as.array(x), np =nuts_params(x), pars =c("alpha_pop[2]","h2_pop[2]","h2_pop[17]","sigma2_r"),off_diag_args =list(size =0.75)) })
$FE
$FS
Comparing estimates
We check the correlations between \(h^2\) and \(I\) estimated population-specific models or with the full model including all populations.
Code
list_df <-lapply(c("I","h2"), function(x){# Extract the median estimates of h2/I in the population-specific models# ======================================================================df_pop_specific_models <-readRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels_FilteredDataset.rds"))) %>%lapply(function(site_i){ site_i %>%lapply(function(pop_i){ broom.mixed::tidyMCMC(pop_i,pars=(x),droppars =NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) }) %>%bind_rows(.id ="pop") }) %>%bind_rows(.id ="site") %>%mutate(model="pop_specific_models") %>%mutate(model_site=paste0(model,"_",site)) %>% dplyr::select(-term)# Extract the median estimates of h2/I in the population-specific models# ======================================================================df <-readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds"))) %>%lapply(function(site_i){ broom.mixed::tidyMCMC(site_i,pars=(paste0(x,"_pop")),droppars =NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) }) %>%bind_rows(.id ="site") %>%mutate(model="full_model") %>%mutate(model_site=paste0(model,"_",site)) %>%mutate(pop=rep(unique(subdata$pop) %>% sort,2)) %>% dplyr::select(-term) %>%bind_rows(df_pop_specific_models)}) %>%setNames(c("I","h2"))list_df[["h2"]] %>%ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site,shape=model)) +geom_pointinterval(position =position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +xlab("") +ylab("Narrow-sense heritability") +theme(axis.text =element_text(size=18),panel.grid.minor.x=element_blank(),panel.grid.major.x=element_blank()) +guides(color=guide_legend(ncol=1))
Archambeau, Juliette, Marta Benito Garzón, Marina de Miguel, Benjamin Brachi, Frédéric Barraquand, and Santiago C González-Martı́nez. 2023. “Reduced Within-Population Quantitative Genetic Variation Is Associated with Climate Harshness in Maritime Pine.”Heredity, 1–11.
Beaton, Joan, Annika Perry, Joan Cottrell, Glenn Iason, Jenni Stockan, and Stephen Cavers. 2022. “Phenotypic Trait Variation in a Long-Term Multisite Common Garden Experiment of Scots Pine in Scotland.”Scientific Data 9 (1): 671.